home *** CD-ROM | disk | FTP | other *** search
/ The X-Philes (2nd Revision) / The X-Philes Number 1 (1995).iso / xphiles / hp48_1 / gausselm < prev    next >
Internet Message Format  |  1995-03-31  |  9KB
  1. Path: seq!spell
  2. From: Robert Brunner <brunner@uirvld.csl.uiuc.edu>
  3. Subject:  v02i018:  gausselm - Gaussian elimination and QR factorization v1.0, Part01/01
  4. Newsgroups: comp.sources.hp48
  5. Followup-To: comp.sys.hp48
  6. Approved: spell@seq.uncwil.edu
  7.  
  8. Checksum:  925599671 (verify with brik -cv)
  9. Submitted-by: Robert Brunner <brunner@uirvld.csl.uiuc.edu>
  10. Posting-number: Volume 2, Issue 18
  11. Archive-name: gausselm/part01
  12.  
  13. BEGIN_DOC gausselm.doc
  14.  
  15. Gaussian elimination and QR factorization
  16. Robert Brunner
  17. brunner@uirvld.csl.uiuc.edu
  18.  
  19. The following directory contains two useful linear algebra routines and
  20. three subroutine functions.
  21.  
  22. QR:  QR decomposes a matrix on the stack, A, into an orthogonal
  23. matrix Q and a upper triangular matrix R such that A=QR, using
  24. Gram-Schmidt orthogonalization.  The program only works for full-rank
  25. matricies.
  26.  
  27. LU:  LU decomposes the matrix from the stack, A, into a permutation
  28. matrix P, a lower triangular matrix L, and an upper triangular matrix
  29. U, such that PA=LU.   Pivoting is only performed when necessary, and
  30. the first non-zero pivot is selected, rather than always choosing the
  31. largest pivot.  By doing this, pivoting is avoided when not needed,
  32. which is useful for solving small problems given for homework in linear
  33. algebra classes.  The method does exhibit poor numerical stability,
  34. so the program is not practical for large matrices.
  35.  
  36. DCMP:  DCMP decomposes a matrix [[row 1][row 2]...[row n]] into a
  37. list of vectors {[row 1][row 2]...[row n]}.  Doing this speeds up
  38. the other two routines.
  39.  
  40. RCMP:  RCMP reverses DCMP.
  41.  
  42. SWRW:  SWRW switches the two specified rows of a matrix decomposed
  43. by DCMP.
  44.  
  45. 'GAUSSELM' BYTES gives checksum #33398d, 1458 bytes
  46.  
  47. Drop me a line if you find this useful.
  48. brunner@uirvld.csl.uiuc.edu
  49.  
  50. END_DOC
  51.  
  52. BEGIN_RPL gausselm.rpl
  53. %%HP: T(3)A(R)F(.);
  54. DIR
  55.   QR
  56.     \<< TRN CONJ DCMP
  57. OVER DUP 2 \->LIST 0
  58. CON { } \-> a n m r q
  59.       \<< a 1 GET DUP
  60. DUP ABS / DUP ROT
  61. DOT 'r' { 1 1 } ROT
  62. PUT 'q' SWAP STO+ 2
  63. n
  64.         FOR i a i
  65. GET DUP \-> ai
  66.           \<< 1 i 1 -
  67.             FOR j
  68. DUP q j GET DUP ai
  69. DOT 'r' { j i } ROT
  70. PUT DUP ROT DOT * -
  71.             NEXT
  72. DUP ABS / DUP ai
  73. DOT 'r' { i i } ROT
  74. PUT 'q' SWAP STO+
  75.           \>>
  76.         NEXT q RCMP
  77. TRN CONJ r
  78.       \>>
  79.     \>>
  80.   LU
  81.     \<< DCMP OVER IDN
  82. DUP 0 CON DCMP
  83. DROP2 SWAP DCMP
  84. DROP2 1 DUP DUP
  85. RCLF \-> a m n l p pr
  86. pc sr flg
  87.       \<<
  88.         WHILE pr m
  89. < pc n \<= AND
  90.         REPEAT 1 CF
  91.           WHILE 1
  92. FC? pc n \<= AND
  93.           REPEAT pr
  94. 'sr' STO
  95.             WHILE 1
  96. FC? sr m \<= AND
  97.             REPEAT
  98.               IF a
  99. sr GET pc GET
  100.               THEN
  101. 1 SF
  102.               ELSE
  103. 'sr' INCR DROP
  104.               END
  105.             END
  106.             IF 1
  107. FC?
  108.             THEN
  109. 'pc' INCR DROP
  110.             END
  111.           END
  112.           IF 1 FS?
  113.           THEN
  114.             IF pr
  115. sr \=/
  116.             THEN
  117. 'p' pr sr SWRW 'l'
  118. pr sr SWRW 'a' pr
  119. sr SWRW
  120.             END a
  121. pr GET DUP pc GET \->
  122. pivr piv
  123.             \<< pr 1
  124. + m
  125.               FOR i
  126. a i GET DUP pc GET
  127. piv / DUP l i GET
  128. pr ROT PUT 'l' i
  129. ROT PUT pivr * -
  130. 'a' i ROT PUT
  131.               NEXT
  132.             \>>
  133.           END 'pr'
  134. INCR 'pc' INCR
  135. DROP2
  136.         END p RCMP
  137. l RCMP m IDN + a
  138. RCMP flg STOF
  139.       \>>
  140.     \>>
  141.   DCMP
  142.     \<< OBJ\-> OBJ\->
  143. DROP DUP2 * OVER 1
  144. - \-> m n mn n1
  145.       \<< 1 m
  146.         FOR i n
  147. \->ARRY mn n1 i * -
  148. ROLLD
  149.         NEXT m
  150. \->LIST m n
  151.       \>>
  152.     \>>
  153.   RCMP
  154.     \<< OBJ\-> OVER
  155. SIZE OBJ\-> DROP OVER
  156. \-> m n st
  157.       \<< 1 m
  158.         START OBJ\->
  159. DROP 'st' n 1 -
  160. STO+ st ROLL
  161.         NEXT 1 n 1
  162. -
  163.         START st
  164. ROLL
  165.         NEXT { m n
  166. } \->ARRY
  167.       \>>
  168.     \>>
  169.   SWRW
  170.     \<< \-> i j
  171.       \<< DUP i GET
  172. SWAP DUP DUP j GET
  173. i SWAP PUT j ROT
  174. PUT
  175.       \>>
  176.     \>>
  177. END
  178. END_RPL
  179.  
  180.  
  181. BEGIN_ASC gausselm.asc
  182. %%HP: T(3)A(D)F(.);
  183. "69A20FF7C780000000403575257540D9D20E16321C432D6E201096D6E2010A6E
  184. 163278BF1D6E2010966C7D1DBBF178BF178BF1D6E2010A66C7D1D6E201096DBB
  185. F1704D1D6E2010A6E0CF1704D1EF53293632B213079000402534D40540D9D20E
  186. 1632B7FC192CF18B9C1B7FC18DBF192CF11C432D6E2010D6D6E2010E6D6E2020
  187. 3747E16329C2A2D6E2010D630132B7FC18DBF145632D6E2020374797632D6E20
  188. 10E69C2A290DA1B4402D6E202037475BCF1C42329C2A2D6E2010E69C2A290DA1
  189. 30132D6E202037475BCF1C423247A20D6E2010D6D6E2010E6B2130900D1EF532
  190. 93632B2130C1100404434D40540D9D20E1632B7FC1B7FC18DBF12ABF1EEDA192
  191. CF19C2A290DA11C432D6E2010D6D6E2010E6D6E2020D6E6D6E2020E613E16329
  192. C2A2D6E2010D60A132D6E201096D6E2010E6900D1D6E2020D6E6D6E2020E613D
  193. 6E201096EEDA190DA10DCF1C4232D6E2010D6387C1D6E2010D6D6E2010E6EF53
  194. 293632B2130CF00020C45520D9D20E163284E20404434D40592CF1CD2D178BF1
  195. 4B2A2681D184E20404434D4053FBF1DBBF184E20404434D4053FBF19C2A278BF
  196. 178BF1916C11C432D6E201016D6E2010D6D6E2010E6D6E2010C6D6E201007D6E
  197. 20200727D6E20200736D6E20203727D6E203066C676E163233032D6E20200727
  198. D6E2010D6EBBE1D6E20200736D6E2010E6CFCE1387E1D5032D9D209C2A25D2C1
  199. 330329C2A2063C1D6E20200736D6E2010E6CFCE1387E1D5032D9D20D6E202007
  200. 2745632D6E2020372797632DCC02330329C2A2063C1D6E20203727D6E2010D6C
  201. FCE1387E1D5032D9D203CE22D6E201016D6E202037276C7D1D6E202007366C7D
  202. 1AFE22D9D209C2A2472C1B21305BF22D9D2045632D6E20203727976324F8028D
  203. BF1B21305DF22B2130496323CE229C2A2063C1AFE22D9D2045632D6E20200736
  204. 976324F8028DBF1B21305DF22B2130496323CE229C2A2313C1AFE22D9D203CE2
  205. 2D6E20200727D6E20203727D9AE1AFE22D9D2045632D6E20100797632D6E2020
  206. 0727D6E2020372784E20403575257545632D6E2010C697632D6E20200727D6E2
  207. 020372784E20403575257545632D6E20101697632D6E20200727D6E202037278
  208. 4E204035752575B21305DF22D6E201016D6E202007276C7D178BF1D6E2020073
  209. 66C7D11C432D6E204007966727D6E2030079667E1632D6E202007279C2A276BA
  210. 1D6E2010D60A132D6E201096D6E201016D6E2010966C7D178BF1D6E202007366
  211. C7D1D6E203007966750FA178BF1D6E2010C6D6E2010966C7D1D6E20200727E0C
  212. F1704D145632D6E2010C697632D6E201096E0CF1704D1D6E204007966727EEDA
  213. 190DA145632D6E20101697632D6E201096E0CF1704D1C4232EF532B21305DF22
  214. 45632D6E20200727976324F80245632D6E20200736976324F8023FBF1B213049
  215. 632D6E20100784E20402534D405D6E2010C684E20402534D405D6E2010D6CD2D
  216. 176BA1D6E20101684E20402534D405D6E203066C676F76C1EF53293632B2130F
  217. A50020152520D9D20E1632293D1E6AA184E20404434D40592CF178BF1ED2A238
  218. 7C14B2A2681D147A20B21301C432D6E201016D6E2010E6D6E2010D6D6E201027
  219. D6E201017E1632D6E2010169C2A26C7D178BF178BF1F1AA150FA178BF1E0CF1E
  220. FFB145632D6E2010279763247A209C2A29C2A2B2130E0CF1704D145632D6E201
  221. 01797632DBBF1B4402ED2A2D6E2010E60A132D6E201096D6E201016D6E201096
  222. 6C7D178BF11C432D6E20201696E16329C2A2D6E2010969C2A290DA10A132D6E2
  223. 010A678BF1D6E201017D6E2010A66C7D178BF1D6E20201696EFFB145632D6E20
  224. 10279763247A20D6E2010A6D6E201096B2130E0CF1704D178BF1E0CF1EFFB1EE
  225. DA190DA1C423278BF1F1AA150FA178BF1D6E20201696EFFB145632D6E2010279
  226. 763247A20D6E201096D6E201096B2130E0CF1704D145632D6E20101797632DBB
  227. F1B4402EF532C4232D6E20101784E20402534D405293D1E6AA1D6E201027EF53
  228. 293632B21306728"
  229. END_ASC
  230.  
  231.  
  232. BYTES: #8276h 1462
  233.  
  234. BEGIN_UU gausselm.uue
  235. begin 644 gausselm
  236. M2%!(4#0X+466*O!_?`@````$4U=25P2=+>!A(\$TTN8"`6EM+A"@YF$CA_O11
  237. MY@(!:<;7T;L?A_MQN!]M+A"@9GP=;2X0D-:['P?4T>8"`6H._'%`'?XUDF,C*
  238. M*S%P"0`$4D--4`2=+>!A(WO/D<(?N,FQ]QS8^Y'"'\$TTN8"`6UM+A#@UN8"1
  239. M`G-T'C:2+"IM+A#0-A`C>\^!O1]4-M+F`@)S='DVTN8"`6[)HI+0&DL$TN8"X
  240. M`G-TM?S!)"/)HM+F`@%NR:*2T!H#,=+F`@)S=+7\P20C="K0Y@(!;6TN$."VK
  241. M$@,)T.%?(SDVLA(#'`%`0#34!$70V0(>-K+W''O/@;T?HOOAWAHI_)$L*@FM#
  242. M$4PC;2X0T-;F`@%N;2X@T.;6Y@(";C$>-I(L*FTN$-`&&B-M+A"0UN8"`6X)W
  243. MT-'F`@)M;FTN(.`6T^8"`6GNK9'0&M#\P20C;2X0T#9X'&TN$-#6Y@(!;OXU[
  244. MDF,C*S'`#P`"3%4"G2W@82-(+D!`--0$E<(?W-)QN!^THF(8'4@N0$`TU`0UM
  245. MOQ^]^X'D`@1$0TU0\_N1+"J'^W&X'QG&$4PC;2X0$-;F`@%M;2X0X-;F`@%L:
  246. M;2X0`-?F`@)P<FTN(``WUN8"`G-R;2XP8,9VYF$C,S#2Y@("<')M+A#0YKL>Q
  247. M;2X@`#?6Y@(!;OSL,7@>73#2V0+)HE(M'#,PDBPJ8,/1Y@("<&-M+A#@QL\>S
  248. M@^?1!2.=+=#F`@)P<E0VTN8"`G-R>3;2S"`S,)(L*F##T>8"`G-R;2X0T,;/5
  249. M'H/GT04CG2TP[")M+A`0UN8"`G-RQM?1Y@("<&/&UZ'O(ITMD"PJ=,*Q$@.UY
  250. M+]+9`E0VTN8"`G-R>39"CR#8^[$2`]4OLA(#E#8R["+)H@(V'/HNTMD"5#;2=
  251. MY@("<&-Y-D*/(-C[L1(#U2^R$@.4-C+L(LFB,C$<^B[2V0+#+M+F`@)P<FTNF
  252. M(#`GUZD>^B[2V0)4-M+F`@%P>3;2Y@("<')M+B`P)X?D`@135U)75#;2Y@(!#
  253. M;'DVTN8"`G!R;2X@,">'Y`($4U=25U0VTN8"`6%Y-M+F`@)P<FTN(#`GA^0"R
  254. M!%-74E<K,5#](FTN$!#6Y@("<'+&UW&X'VTN(``W9GP=P332Y@($<&EV<FTN3
  255. M,`"79N=A(VTN(``GERPJ9ZO1Y@(!;:`QTN8"`6EM+A`0UN8"`6G&UW&X'VTNN
  256. M(``W9GP=;2XP`)=F5_`:A_O1Y@(!;&TN$)!F?!UM+B``)^?`'P?4064C;2X0E
  257. MP)9G(VTN$)#FP!\'U-'F`@1P:79R[JV1T!I4-M+F`@%A>3;2Y@(!:0[\<4`=,
  258. M3#+B7R,K,5#](E0VTN8"`G!R>39"CR!4-M+F`@)P8WDV0H\@\_NQ$@.4-M+FO
  259. M`@%P2"Y`(#74!-7F`@%L2"Y`(#74!-7F`@%MW-)QMAIM+A`0AN0"!%)#35!M[
  260. M+C!@QG;V9QS^-9)C(RLQ\%H``E%2`ITMX&$CDM/AIAI(+D!`--0$E<(?A_OAB
  261. M+2J#QT$K*H;10:<"*S$03"-M+A`0UN8"`6YM+A#0UN8"`7)M+A`0YV$C;2X0I
  262. M$)8L*L;7<;@?A_OQH1H%KW&X'P[\X?\;5#;2Y@(!<GDV0J<"R:*2+"HK,>#`B
  263. M'P?4064C;2X0$)=G([W[L40@WJ+2Y@(!;J`QTN8"`6EM+A`0UN8"`6G&UW&X2
  264. M'\$TTN8"`F%I'C:2+"IM+A"0EBPJ":T!&B-M+A"@=K@?;2X0$-?F`@%JQM=QQ
  265. MN!]M+B`0EN;_&U0VTN8"`7)Y-D*G`FTN$*#6Y@(!:2LQX,`?!]1QN!\._.'_$
  266. M&^ZMD=`:3#)RN!\?JE'P&H?[T>8"`F%I_K]!92-M+A`@EV<C="K0Y@(!:6TNP
  267. M$)"V$@,._'%`'50VTN8"`7%Y-M*['TL$XE\C3#+2Y@(!<4@N0"`UU`0E.1UNM
  268. .JM'F`@%R_C628R,K,0`KM
  269. ``
  270. end
  271. END_UU
  272.